PROGRAM pondSim
    IMPLICIT NONE
    !PARAMETERS
    DOUBLEPRECISION::tol1=5e-3 !max change in ft during single timestep (dictates max delta.t)
    DOUBLEPRECISION::tol2=1e-3 !min change in ft during single timestep (dictates condition where delta.t increased)
    DOUBLEPRECISION::tol3=1 !dictates min delta.t
    DOUBLEPRECISION::tFactor=1.2 !factor by which delta.t increased/decreased
    DOUBLEPRECISION::mgdToCfs=1.55 !conversion factor
    DOUBLEPRECISION::pi=3.14
    DOUBLEPRECISION::grav=32.2 !ft/s^2

    !VARIABLES USED IN SUBROUTINE readInput
    INTEGER::i,j !counter variables
    INTEGER::nDays !dys
    INTEGER::nPonds !number of ponds
    INTEGER::nCons !number of connections
    INTEGER::dummy !just a placeholder for unused input info
    INTEGER::nMaxFlows,nMaxFlowCons !number of max flow specs, # cons associated w/ given flow
    INTEGER,DIMENSION(3)::date,starttime,stoptime !system date and time vectors
    INTEGER,DIMENSION(:),ALLOCATABLE::setPond !0=none, 1=specify height, 2=specify flow
    INTEGER,DIMENSION(:),ALLOCATABLE::setCon !0=none, 1=specify flow
    INTEGER,DIMENSION(:),ALLOCATABLE::receivesPlantInflow !0=no, 1=yes
    INTEGER,DIMENSION(:),ALLOCATABLE::upstreamPond !upstream pond number
    INTEGER,DIMENSION(:),ALLOCATABLE::dnstreamPond !dnstream pond number
    INTEGER,DIMENSION(:),ALLOCATABLE::conType !1=pipe, 2=rect weir, 3=v-weir, 4=pump
    INTEGER,DIMENSION(:),ALLOCATABLE::conAdjust !0=no, 1=yes
    INTEGER,DIMENSION(:,:),ALLOCATABLE::maxFlowCons !columns = cons associated w/ max flow
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::hTemp !mgd
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::maxFlow !mgd
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::pondMinFlow !mgd
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::pondMaxFlow !mgd
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::plantInflow !mgd
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::precip,evap !inches per day
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::pondArea !ft^2
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::initialHt,minHt,maxHt !ft above msl
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::pondTargetValue !ft above msl for height, !mgd for flow
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::upPipeLength,dnPipeLength !ft
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::upPipeFric,dnPipeFric !0-1
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::upPipeDiam,dnPipeDiam !ft
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::upPipeSluice,dnPipeSluice !0-100 (percent open)
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::weirWidth,weirHeight !ft,ft above msl
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::weirFric !friction factor 0-1
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::weirOverflow !ft above msl
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::conTargetValue !mgd for flow
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::setConMinHt !mgd for flow
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::pumpCapacity !mgd for flow
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::pumpHoursOn !total on-time (hrs)
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::pumpLagTime !hrs before on-time
    CHARACTER(6),DIMENSION(:),ALLOCATABLE::pondNames
    CHARACTER(6),DIMENSION(:),ALLOCATABLE::conNames
    CHARACTER(100),DIMENSION(:),ALLOCATABLE::maxFlowNames

    !VARIABLES USED IN SUBROUTINE runSim
    INTEGER::step,day,hour,previousHour,adjustInterval !counters, adjustInterval=hrs between weir adjustments
    INTEGER::maxSteps !max # steps = total seconds in simulation
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::totalTime !net simulated time (secs) associated w/ each timeStep
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::pondFlow !mgd
    DOUBLEPRECISION,DIMENSION(:),ALLOCATABLE::conFlow !mgd
    DOUBLEPRECISION,DIMENSION(:,:),ALLOCATABLE::h !pond levels in ft above msl
    DOUBLEPRECISION,DIMENSION(:,:),ALLOCATABLE::K !RK dummy variables
    DOUBLEPRECISION::deltaT !timestep length   
    CHARACTER(100)::fmt

    CALL readInput !READS INPUT FILE
    CALL runSim !RUN SIMULATION
    CALL printOutput !WRITE OUTPUT FILE

    CONTAINS
    
    SUBROUTINE readInput
        IMPLICIT NONE
        
        !NOW READ VARIABLES FROM FILE
        OPEN(111,FILE='pondSim.input',STATUS='OLD')
        
        !READ DAILY DATA
        READ(111,*) nDays
        ALLOCATE(precip(nDays),evap(nDays),plantInflow(nDays))
        DO i=1,nDays
            READ(111,*) plantInflow(i),precip(i),evap(i)
        END DO
        
        !READ POND PARAMETERS
        READ(111,*) nPonds
        ALLOCATE(setPond(nPonds),initialHt(nPonds),minHt(nPonds),maxHt(nPonds),hTemp(nPonds),&
        pondTargetValue(nPonds),pondArea(nPonds),pondNames(nPonds),receivesPlantInflow(nPonds),&
        pondMinFlow(nPonds),pondMaxFlow(nPonds))
        DO i=1,nPonds !ONE LOOP PER POND
            READ(111,*) dummy
            READ(111,*) pondNames(i)
            READ(111,*) pondArea(i)
            READ(111,*) receivesPlantInflow(i)
            READ(111,*) minHt(i)
            READ(111,*) maxHt(i)
            READ(111,*) initialHt(i)
            READ(111,*) pondMinFlow(i)
            READ(111,*) pondMaxFlow(i)
            READ(111,*) setPond(i)
            IF (setPond(i)==1 .OR. setPond(i)==2) READ(111,*) pondTargetValue(i)
        END DO

        !READ CONNECTION PARAMETERS
        READ(111,*) nCons
        ALLOCATE(setCon(nCons),setConMinHt(nCons),upstreamPond(nCons),dnstreamPond(nCons),conType(nCons),&
                upPipeLength(nCons),dnPipeLength(nCons),upPipeFric(nCons),dnPipeFric(nCons),&
                upPipeDiam(nCons),dnPipeDiam(nCons),upPipeSluice(nCons),dnPipeSluice(nCons),&
                weirWidth(nCons),weirHeight(nCons),weirFric(nCons),weirOverflow(nCons),conTargetValue(nCons),&
                pumpCapacity(nCons),pumpHoursOn(nCons),pumpLagTime(nCons),conNames(nCons),conAdjust(nCons))
        !INITIALIZE A FEW PARAMETERS
        conAdjust=0
        setCon=0
        setConMinHt=0
        conAdjust=0

        DO i=1,nCons !ONE LOOP PER CONNECTION
            READ(111,*)dummy
            READ(111,*) conNames(i)
            READ(111,*) upstreamPond(i)
            READ(111,*) dnstreamPond(i)
            READ(111,*) conType(i)
            !READ PIPE INFO
            IF (conType(i)==1) THEN
                READ(111,*) upPipeLength(i)
                READ(111,*) upPipeDiam(i)
                READ(111,*) upPipeFric(i)
            END IF
            !READ WEIR INFO
            IF (conType(i)==2 .OR. conType(i)==3) THEN
                READ(111,*) upPipeLength(i)
                READ(111,*) upPipeDiam(i)
                READ(111,*) upPipeFric(i)
                READ(111,*) upPipeSluice(i)
                READ(111,*) dnPipeLength(i)
                READ(111,*) dnPipeDiam(i)
                READ(111,*) dnPipeFric(i)
                READ(111,*) dnPipeSluice(i)   
                READ(111,*) weirWidth(i) 
                READ(111,*) weirHeight(i)
                READ(111,*) weirFric(i)
                READ(111,*) weirOverflow(i)
                READ(111,*) conAdjust(i)
                READ(111,*) setCon(i)
                IF (setCon(i)==1) THEN
                    READ(111,*) conTargetValue(i)
                    READ(111,*) setConMinHt(i)
                END IF
            END IF
            !READ PUMP INFO
            IF(conType(i)==4) THEN
                READ(111,*) pumpCapacity(i)
                READ(111,*) pumpHoursOn(i)
                READ(111,*) pumpLagTime(i)
        END IF
        END DO
        
        !READ MAX FLOW INFO
        READ(111,*) nMaxFlows
        ALLOCATE(maxFlow(nMaxFlows),maxFlowNames(nMaxFlows),maxFlowCons(nMaxFlows,nCons))
        maxFlowCons=0 !INITIALIZE TO ZERO
        DO i=1,nMaxFlows
            READ(111,*)maxFlowNames(i)
            READ(111,*)maxFlow(i)
            READ(111,*)nMaxFlowCons
            DO j=1,nMaxFlowCons
                READ(111,*) dummy
                maxFlowCons(i,dummy)=1
            END DO
        END DO
            
        CLOSE(111)
    END SUBROUTINE readInput   
        
    SUBROUTINE runSim
        IMPLICIT NONE
   
        !GET SYSTEM DATE AND TIME
        CALL IDATE(date)
        CALL ITIME(starttime)
 
        !ALLOCATE ARRAYS
        maxSteps=nDays*24*60*60
        ALLOCATE(totalTime(maxSteps),h(nPonds,maxSteps),K(nPonds,3))
        ALLOCATE(pondFlow(nPonds),conFlow(nCons))
        
        !OPEN LOG/OUTPUT FILES FOR WRITING
        OPEN(111,FILE='pondSim.output') !output filename
        OPEN(112,FILE='pondSim.log') !logfile name
        OPEN(113,FILE='pondSim.diagnostics') !logfile name
        fmt='(A,I2,A,I2,A,I4,A,I2,A,I2,A,I2)'
        WRITE(112,fmt)&
            'Begin Simulation. Date: ',date(2),'/',date(1),'/',date(3),&
            '  Hour: ',starttime(1),':',starttime(2),':',starttime(3)
        
        !INITIALIZE COUNTERS/VARIABLES
        step=1
        totalTime(1)=0d0
        adjustInterval=6 !hrs
        previousHour=0
        deltaT=1d0 !sec
        h(:,1)=initialHt
        
        !BEGIN SIMULATION
        DO 
            step=step+1
            totalTime(step)=totalTime(step-1)+deltaT
            day=1+INT(totalTime(step)/24d0/60d0/60d0)
            hour=1+INT(totalTime(step)/60d0/60d0)
            
            !CHECK FOR EXIT CRITERIA
            IF (totalTime(step)>=maxSteps) THEN
                CALL ITIME(stoptime)
                fmt='(A,I2,A,I2,A,I2)'
                WRITE(112,fmt)&
                    'Simulation Succesfull! Hour: ',stoptime(1),':',stoptime(2),':',stoptime(3)
                RETURN
            END IF

            !MAKE ADJUSTMENTS AS NECCESSARY
            IF ((hour>previousHour .AND. MOD(hour,adjustInterval)==0) .OR. step<=2) THEN
                CALL adjustCons
                deltaT=1 !RESET TIMESTEP TO 1 SECOND
            END IF
            previousHour=hour

            !RK STEP #1
            CALL findFlow(h(:,step-1))   
            pondFlow=pondFlow+precip(day)-evap(day)
            K(:,1)=deltaT*(pondFlow/pondArea)
            hTemp=h(:,step-1)+(1d0/3d0)*K(:,1)

            !RK STEP #2
            CALL findFlow(hTemp)
            pondFlow=pondFlow+precip(day)-evap(day)
            K(:,2)=deltaT*(pondFlow/pondArea)
            hTemp=h(:,step-1)+(2d0/3d0)*K(:,2)

            !RK STEP #3
            CALL findFlow(hTemp)
            pondFlow=pondFlow+precip(day)-evap(day)
            K(:,3)=deltaT*(pondFlow/pondArea)
            hTemp=h(:,step-1)+(1d0/4d0)*K(:,1)+(3d0/4d0)*K(:,3)
    
            !STORE HEADS AND UPDATE TOTAL TIME
            h(:,step)=hTemp
            totalTime(step)=totalTime(step-1)+deltaT
 
            !IMPLEMENT VARIABLE TIMESTEP
            IF ((MOD(totalTime(step)/60d0/60d0,DBLE(adjustInterval))>(10d0/60d0)) &
            .AND. ANY(ABS(hTemp-h(:,step-1)) >= tol1) ) THEN
                !CHANGE IN HEADS TO LARGE, REDUCE DELTA-T
                IF ((deltaT/tFactor)>tol3) THEN
                    deltaT=deltaT/tFactor
                    step=step-1
                    CYCLE
                ELSE
                    WRITE(111,*) '0' ! #SIMULATION UNSUCCESSFUL'
                    WRITE(112,*) 'ERROR: min timestep not small enough to capture change in levels'
                    WRITE(112,*) 'step: ',step
                    WRITE(112,*) 'hour: ',totalTime(step)/60d0/60d0
                    WRITE(112,*) 'timestep: ',deltaT,' tol: ',tol1
                    DO i=1,nPonds
                        WRITE(112,*) pondNames(i),h(i,step-1)-hTemp(i)
                    END DO
                    STOP
                END IF
            ELSE IF ( ALL(ABS(hTemp-h(:,step-1)) <= tol2)  ) THEN
                !CHANGE IN HEADS SMALL ENOUGH TO INCREASE DELTA-T
                deltaT=deltaT*tFactor
            END IF

            !IF TIMESTEP OK, CHECK FOR LOW/HIGH LEVELS
            IF (ANY(h(:,step)<minHt)) THEN
                WRITE(111,*) '0' ! #SIMULATION UNSUCCESSFUL'
                DO j=1,nPonds
                    IF (h(j,step)<minHt(j)) THEN
                        WRITE(112,*) 'ERROR: low level in pond: ',pondNames(j)
                        WRITE(112,*) 'minlevel: ',minHt(j),'level: ',h(j,step)
                        WRITE(113,'(42F20.10)') totalTime(step),h(:,step),conFlow 
                    END iF
                END DO
                STOP
            ELSE IF (ANY(h(:,step)>maxHt)) THEN
                WRITE(111,*) '0' ! #SIMULATION UNSUCCESSFUL'
                DO j=1,nPonds
                    IF (h(j,step)>maxHt(j)) THEN
                        WRITE(112,*) 'ERROR: high level in pond: ',pondNames(j)
                        WRITE(112,*) 'maxlevel: ',maxHt(j),'level: ',h(j,step) 
                        WRITE(113,'(42F20.10)') totalTime(step),h(:,step),conFlow
                    END iF
                END DO
                STOP
            END IF

            !CHECK FOR HIGH FLOWS
            DO i=1,nMaxFlows
                dummy=0
                DO j=1,nCons
                    IF(maxFlowCons(i,j)==1) dummy=dummy+conFlow(j)
                END DO    
                IF (dummy>maxFlow(i)) THEN
                    WRITE(111,*) '0'! #SIMULATION UNSUCCESSFUL'
                    WRITE(112,*) 'ERROR: High Flow Through: ',maxFlowNames(i)
                    STOP
                END IF
            END DO      
            
            WRITE(113,'(42F20.10)')totalTime(step),h(:,step),conFlow
        END DO
        
        CLOSE(111)
        CLOSE(112)
        CLOSE(113)
    END SUBROUTINE runSim

    SUBROUTINE adjustCons
        IMPLICIT NONE
        DOUBLEPRECISION::dummydble       
 
        !ADJUSTS 'AUTO-WEIRS'
        DO i=1,nCons
            IF ( setCon(i)==1 .AND. h(upstreamPond(i),step-1)>setConMinHt(i) ) THEN 
                IF (conType(i)==2) THEN
                    weirHeight(i) = h(upstreamPond(i),step-1)-((3d0/10d0)*conTargetValue(i)*mgdToCfs/weirWidth(i))**(2d0/3d0)
                ELSE IF (conType(i)==3) THEN
                    weirHeight(i) = h(upstreamPond(i),step-1)-((2d0/5d0)*conTargetValue(i)*mgdToCfs)**(2d0/5d0)
                END IF
            END IF
        END DO

        !NEXT LOOP ADJUSTS 'AUTO-PONDS'
        !PRECEDENCE:    can always lower head, but warn if resulting dnstream flow is large
        !               can always raise head, but warn if resulting dnstream flow is small
        !               cannot increase flow rate if resulting upstream head is too low
        !               cannot decrease flow rate if resulting upstream head is too high
        DO i=1,nPonds

            !IF POND LEVEL IS SPECIFIED, ADJUST ALL OUTGOING WEIRS
            IF(setPond(i)==1) THEN
                !for each set of weirs (set = leading to a common downstream pond)
                DO j=1,nCons
                    dummydble=weirHeight(j)
                    IF (upstreamPond(j)==i .AND. conAdjust(j)==1) THEN
                        weirHeight(j)=pondTargetValue(i)
                        
                        !RE-ADJUST WEIR IF BELOW PRESCRIBED MIN
                        IF ( weirHeight(j) < setConMinHt(j) ) weirHeight(j)=setConMinHt(j)
                        IF ( weirHeight(j) < minHt(i)) weirHeight(j)=minHt(i)
 
                        !WRITE WARNING TO LOGFILE
                        IF ((weirHeight(j)-dummydble)<(-1d-3)) THEN
                            WRITE(112,*)'WARNING: Hour: ',totalTime(step-1)/60d0/60d0, ' Raising flow through con: ',conNames(j) 
                        ELSE IF ((weirHeight(j)-dummydble)>(1d-3)) THEN
                            WRITE(112,*)'WARNING: Hour: ',totalTime(step-1)/60d0/60d0, ' Reducing flow through con: ',conNames(j) 
                        END IF 
                    END IF
                END DO

            !IF POND FLOW IS SPECIFIED, ADJUST ALL INCOMING WEIRS
            ELSE IF(setPond(i)==2) THEN
                dummy=COUNT(dnstreamPond==i .AND. conAdjust==1) !dummy = # of adjustable weirs whose dnstream pond = i
                DO j=1,nCons
                    dummydble=weirHeight(j)
                    IF (dnstreamPond(j)==i .AND. conAdjust(j)==1) THEN
                        IF (conType(j)==2) THEN
                            weirHeight(j) = h(upstreamPond(j),step-1)-((3d0/10d0)*pondTargetValue(i)/DBLE(dummy)&
                            *mgdToCfs/weirWidth(j))**(2d0/3d0)
                        ELSEIF (conType(i)==3) THEN
                            weirHeight(j) = h(upstreamPond(j),step-1)-((2d0/5d0)*pondTargetValue(i)/DBLE(dummy)&
                            *mgdToCfs)**(2d0/5d0)
                        END IF
                        
                        !RE-ADJUST WEIR IF BELOW PRESCRIBED MIN
                        IF ( weirHeight(j) < setConMinHt(j) ) weirHeight(j)=setConMinHt(j)
                        IF ( weirHeight(j) < minHt(i)) weirHeight(j)=minHt(i)
                                    
                        !WRITE WARNING TO LOGFILE
                        IF ((weirHeight(j)-dummydble)<(-1d-3)) THEN
                            WRITE(112,*)'WARNING: Hour: ',totalTime(step-1)/60d0/60d0, ' Raising flow through con: ',conNames(j)
                        ELSE IF ((weirHeight(j)-dummydble)>(1d-3)) THEN
                            WRITE(112,*)'WARNING: Hour: ',totalTime(step-1)/60d0/60d0, ' Reducing flow through con: ',conNames(j)
                        END IF
                    END IF
                END DO
            END IF
        END DO
    END SUBROUTINE adjustCons
    
    SUBROUTINE findFlow(hh)
        IMPLICIT NONE
        DOUBLEPRECISION,DIMENSION(:),INTENT(IN)::hh
        DOUBLEPRECISION::pipeArea
        conFlow=0d0 !INITIALIZE FLOWS TO ZERO
        pondFlow=0d0
        DO i=1,nCons
            !CALCULATE FLOW ACROSS CONNECTION
            IF (conType(i)==1) THEN !PIPE
                pipeArea=(pi*(upPipeDiam(i)/2)**2)
                conFlow(i)=pipeArea*SQRT(2*grav*ABS(hh(upstreamPond(i))-hh(dnstreamPond(i))))
                !ADJUST SIGN FOR REVERSE FLOW THROUGH PIPE
                IF (hh(upstreamPond(i))<hh(dnstreamPond(i))) conFlow(i)=(-1)*conFlow(i)
            ELSE IF (conType(i)==2) THEN !RECTANGULAR WEIR
                IF(hh(upstreamPond(i))>weirHeight(i)) THEN
                    conFlow(i)=(10d0/3d0)*weirWidth(i)*(hh(upstreamPond(i))-weirHeight(i))**(1.5d0)
                END IF
            ELSE IF (conType(i)==3) THEN !V-NOTCH WEIR
                IF (hh(upstreamPond(i))>weirHeight(i)) THEN
                    conFlow(i)=(2.5d0)*(hh(upstreamPond(i))-weirHeight(i))**(2.5d0)
                END IF
            ELSE IF (conType(i)==4) THEN !PUMP
                IF (totalTime(i)/60d0/60d0 > pumpLagTime(i)) THEN
                    IF (totalTime(i)/60d0/60d0 < (pumpLagTime(i) + pumpHoursOn(i))) THEN
                        conFlow(i)=pumpCapacity(i)*mgdToCfs
                    END IF
                END IF
            END IF
            !FIND NET PONDFLOW
            pondFlow(upstreamPond(i))=pondFlow(upstreamPond(i))-conFlow(i) !SUBTRACT FOR OUTFLOW
            pondFlow(dnstreamPond(i))=pondFlow(dnstreamPond(i))+conFlow(i) !ADD FOR INFLOW
        END DO

        !ADD PLANT INFLOW TO APPROPRIATE POND, and PRECIP/EVAP TO ALL PONDS
        DO i=1,nPonds
            IF (receivesPlantInflow(i)==1) THEN
                pondFlow(i)=pondFlow(i)+plantInflow(day)*mgdToCfs
            END IF
            pondFlow(i)=pondFlow(i)+(precip(day)-evap(day))*pondArea(i)/12d0/24d0/60d0/60d0 
            !divide by 12*24*60*60: in/day -> ft/sec
        END DO
    END SUBROUTINE findFlow
    
    SUBROUTINE printOutput
        IMPLICIT NONE
        OPEN(111,FILE='pondSim.output') !output file
        WRITE(111,*) '1' ! #SIMULATION SUCCESSFUL' 
        DO i=1,nPonds
            WRITE(111,*) pondNames(i)
            WRITE(111,*) minHt(i)
            WRITE(111,*) maxHt(i)
            WRITE(111,*) initialHt(i)
            WRITE(111,*) MAXVAL(h(i,1:step-1))
            WRITE(111,*) h(i,step-1)
        END DO
        CLOSE(111)
    END SUBROUTINE printOutput

END PROGRAM pondSim
        
          
